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Abstract 

J. Math. Phys., 42(2001), 200-212. The Fokker-Planck equation 
has been very useful for studying dynamic behavior of stochastic differ- 
ential equations driven by Gaussian noises. In this paper, we derive a 
Fractional Fokker-Planck equation for the probability distribution of par- 
ticles whose motion is governed by a nonlinear Langevin-type equation, 
which is driven by a non-Gaussian Levy-stable noise. We obtain in fact 
a more general result for Markovian processes generated by stochastic 
differential equations. 

PACS numbers: 05.40+j,05.60+w, 66.10Cb, 05.70 a. 
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1 Introduction and statement of the prob- 
lem 

The Fokker-Planck equation is one of the most celebrated equations in 
Physics, since it has been very useful for studying dynamic behavior of 
stochastic differential equations driven by Gaussian noises. There has 
recently been a mushrooming interest [1, 2, 3, 4, 5, 6] in the fact that 
the probability density of a (linear) Levy motion satisfies a generalized 
Fokker-Planck equation involving fractional orders of differentiation. This 
essentially corresponds to a re-interpretation of the characteristic function 
of a Levy motion and significant applications seem to require its nonlinear 
generalization. 

We therefore consider the following nonlinear Langevin-like equation 
for a stochastic (real) quantity X(t) : 

dX(t) = m(X(t),t) dt + o(X(t),t) dL (1) 

where the driving source is a Levy stable motion L, i.e. a motion (e.g. 
[7]) whose increments SL are stationary and independent for any time 
lag St and correspond to independent, identically distributed Levy stable 
variables [9, 10, 11, 12, 13]. Let us recall that a Levy stable motion 
is defined, as are its increments, by four parameters: its Levy stability 
index a (0 < a < 2) , its skewness (3 ( — 1 < /3 < 1), its center ^8t 
and its scale parameter D St (D > 0). Brownian motion corresponds 
to the limit case a = 2, which also implies (3 = 0, and to the 'normal' 
diffusion law. The variance Var[X(t) — X(to)] of the distance traveled by 
a brownian particle, is twice its scale parameter and therefore yields the 
classical Einstein relation: Var[X(t) — X(to)] = 2D(t — to). 

The linear case, which is the unique case studied until now, corre- 
sponds to: 

m{X(i),t)=m = Const., a(X(t),i) = a = Const. (2) 

X(t) — X(to) is also a Levy motion which has the same Levy stability 
index a, but with a possible different center or trend (when m =fc 0) and 
scale or amplitude (when a =fc 1). 

In the nonlinear cases, a(X(t),t) > and m(X(t),t) are nonlinear 
functions of X(t) and t, which satisfy certain regularity constraints to be 
discussed later. We claim that: 

Proposition 1 The transition probability density: 

Vt > to : p(x, t\x ,t ) = Pr(X(t) = x\X(t ) = x ) (3) 

corresponding to the nonlinear stochastic differential equation 
(Eq.l), with a / 1 or /3 = 0, is the solution of the following 
Fractional Fokker-Planck equation: . 
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d d 

—p(x,t\x ,t ) = - — (ja(x,t) + m(x,t))p(x,t\x ,t ) 
-D[(-A) a/2 + ^(a)^(-A) {a - 1)/2 ]<T(x,t) a p(x,t\x ,to) (4) 
where co(a) is defined by: 

a ^ 1 : u)(a) = tan^- (5) 

and where the fractional powers of the Laplacian A will be discussed in 
Sect. 5. Proposition 1 and Eq.4 will be established for scalar processes 
(i.e. A = -g^z) and its extension to vector processes will be discussed and 
presented in Sect. 7. 

This Fractional Fokker-Planck equation will be established with the 
help of the much more general proposition: 

Proposition 2 The inverse Fourier transform of the second 

characteristic function or cumulant generating function of the 
increments of a Markov process X(t) generates by convolution 
the Fokker-Planck equation of evolution of its transition prob- 
ability p(x,t\xo, to). 

We will demonstrate this proposition in a straightforward, yet rigorous 
way. More precisely, we will establish the following: 

^(x,t\x ,to) = J dy^(x -y\y,t)p(y,t\x ,t ) (6) 

where K is the inverse Fourier transform of the cumulant generating func- 
tion of the increments. Its arguments will become explicit in Sect. 2. 

This not only holds for processes with stationary and independent 
increments, as in the linear case (Eq.2) but for any Markov process, in- 
cluding those defined by the non-linear Langevin-like equation (Eq.l with 
m Const., a / Const.). As a consequence of Eq.6, we will demonstrate 
the following: 

Proposition 3 The Kramers-Moyal coefficients A n of the Fokker- 

Planck equation of a Markov process X(t): 

^(x,t\x ,t ) = ^ ^^[An(x,t)p(x,t\x ,to)] (7) 

neJ 

are directly related to the cumulants C n of the increments: 

A n {x,t) = i~^-C n {x,t) (8) 
n! 

where the set J of the indices n is {1,2} in the most classical case (e.g. 
[8], which is a particular case of J C N which corresponds to an analytic 
expansion of cumulants. 
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We will demonstrate this property (Prop. 3), which at best is only 
mentioned in few standard text books on the (classical) Fokker-Plank 
equation, as well as its generalization for non analytic cumulant expan- 
sions, i.e. there are non integers indices n £ J. This latter property, 
discussed in Sect. 5, will be exploited in Sect. 6 in order to derive Prop.l 
with J = {l,a},0 < a < 2. 

2 The cumulant generating function of 
the increments 

The first and second (conditional) characteristic functions are respectively 
the moment generating function Zx(k,t — to\xo,to) and the cumulant 
generating function Kx(k,t — to\xo,to), associated with the transition 
probability p(x, t\xo, to) of a process X(t). They are defined by the Fourier 
transform of the latter, with k being the conjugate variable of x — xq : 



F[p(x,t\x ,t )] = Z x (k,t - t \x ,t ) (9) 
= exp(K x (k,t — to\xo,t )) (10) 
ee E[exp(ik(X{t)-X ))\X(to)] = X ] (11) 

where E[-\-] denote the conditional mathematical expectation, F and F^ 1 
respectively the Fourier-transform and its inverse: 

/oc 
dx exp(ikx)f(x) (12) 
-oo 

f°° dk 

F^[f] = f{x) = / — exp(-ikx)f(k) (13) 

J — oc 

The corresponding quantities for increments 5X(5t) = X(t+St)—X(t), 
corresponding to a given time lag St > 0, are defined in a similar way: 

F\p(x + 5x,t + St\x,t)] = 5Z x {k,St\x,t) (14) 
= exp(8K x {k,St\x,t)) (15) 
= E[exp(ik(X{t + 5t) - X))\X(t) = X\16) 

where k is the conjugate variable of 5x. The cumulants of the increments 
C n are the coefficients of the Taylor expansion of SKx ■ 

SK x (k,8t\x,t) = St^2 ^ r C n (x,t)+o(8t) (17) 

As already mentioned, the classical case corresponds to an analytic 
expansion of SKx, i.e. J C N, whereas we will be interested by a non- 
analytic expansion J — {1, a}. 
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3 Processes with stationary and indepen- 
dent increments 



Let us first consider the simple sub-case of a process with stationary and 
independent increments. It corresponds to C n (x,t) = C n — Const, in 
Eqs.7, 17 and as already discussed in Sect. 1, it includes the linear case 
(Eq.2) of the Langevin-like equation (Eq.l). 

However we believe that the following derivation is not only somewhat 
pedagological on the role of the characteristic functions for the nonlinear 
case, but terser than derivations previously presented for the linear case. 

The stationarity of the increments implies that the transition proba- 
bility depends only on the time and space lags, i.e.: 

p(x, t\x , t ) = p(x - x , t - t ) (18) 

and similarly, the characteristic functions of the increments are no longer 
conditioned, for instance: 

Zx(k,t-to\x ,t ) = Z x (k,t-t ) (19) 
Kx(k,t-to\x ,t ) = K x (k,t-t Q ) (20) 

On the other hand, the independence of the increments implies that 
the transition probabilities satisfy a convolution (over any possible inter- 
mediate position y) for any given time lag 8t: 

\/8t > : p(x — xo,t + St — to) — J dy p(x — y, St)p(y — xq, t — to) (21) 

and the corresponding characteristic functions merely factor (resp. add). 
Therefore, we have: 

Z x (k,t + St - t ) - Z x (k,t- t ) = Z x {k, t - t )(SZ x (k, St) - 1) (22) 
This in turn leads to: 

Z x (k,t + St-t ) - Z x (k,t-t ) = Zx(k,t~to)SK x (k,8t) + o(St) (23) 
Its inverse Fourier transform yields: 



p(x,t+St\xo,to)-p(x,t\x ,t ) = J dyF 1 [SK x (k, St)]p(y-xo,t-t Q )+o(St) 

(24) 

This demonstrates (in the limit St — » 0) Prop. 2 and Eq.6, as well as 
Prop. 3, since Eq.24 corresponds, with the help of Eq.17, to: 
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p(x,t+6t\xo,to)-p(x,t\x ,to) = Sty^JCn j— / dy5^ l } y p{y,t\x ,t )]+o{5t) 

(25) 

where 8" denotes the n th derivative of the Dirac function. Therefore, we 
obtain: 

—p(x,t\x ,t ) =} j A n -Q^p(x,t\x ,to) (26) 
which corresponds to the linear case of Eq.7. 

4 More general Markov processes 

In the case of a Markov process which does not have stationary and in- 
dependent increments, there is no longer a simple convolution equation 
(Eq. 21) of the transition probabilities, nor a simple factorization of char- 
acteristic functions (Eq.22). However, the former satisfies a generalized 
convolution equation which corresponds to the Chapman-Kolmorogorov 
identity [12] valid for any Markov process X(t): 



V5t>0:p(x,t + 5t\x ,t )) = J dy p(x,t + 5t\y,t)p(y,t\x ,t ) (27) 

which indeed reduces to a mere convolution (Eq. 21) in the case of pro- 
cesses with stationary and independent increments. This identity can be 
written under the equivalent form: 

p(x,t + 5t\x ,t ) = j dy J ^e- ikv+SKxlkM) p(y,t\x ,t ) (28) 
Noting that we have: 

p(x,t\x ,t ) = J dy p(y,t\x ,to) J ^ e ~ iky (29) 

we obtain: 

p(x,t+St\x ,to)-p(x,t\x ,to) — St J dyF~ 1 [8Kx(k,St\y,t)]p(y,t\x Q ,to)+o(St) 

(30) 

In the limit St — > 0, this corresponds to Prop. 2 and Eq. 6. When 
J C N, it yields with the help of Eq.17: 



dpi r./|.n,.f ) StJ2 I d vfc-v\^f- C "(l > ni>(U-t\->-« t„)]+o(,)f) CM.) 

nSJV J 

The limit St — > corresponds to Eq.7 and demonstrates Prop. 3 for 
any (classical) Markow process. 
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5 Extension to fractional orders 



In the two previous sections (Sects. 3- 4), the fact that the indices n G J 
should be integers intervene at best only in the correspondence between 
(integer order) differentiation J^j- (in Eq. 7) and powers of the conjugate 
variable fe™ (in Eq. 17). However, by the very definition of fractional 
differentiation (e.g. [14]), this correspondence holds also for non integer 
orders. However, there is not a unique definition of fractional differentia- 
tion and therefore, as discussed in some details in [6]), we cannot expect 
to have a unique expression of the Fractional Fokker-Planck equation. 

Since it will be sufficient for the following to consider an expansion 
of the characteristic function involving fractional powers of only the wave 
number |fc|, it is interesting to consider Riesz's definition of a fractional 
differentiation. Indeed, the latter corresponds to consider fractional pow- 
ers of the Laplacian: 

-(-A) a/2 f(x) = F- 1 [\k\ a f(k)] (32) 

which has furthermore the advantage of being valid for the vector cases. 
However, we will see in Sect. 7 that in general it does not apply in a 
straightforward manner for d-dimensional stable Levy motions. Indeed 
the latter introduces rather (one-dimensional) directional Laplacians, i.e. 
(one-dimensional) Laplacians along a given direction u (| u \= 1) : 

-(-A^ 2 f(x) = F- 1 [|(fc,u)| a /(fc)] (33) 
where(.,.) denotes the scalar product. On the other hand, it will be useful 
to consider the fractional power of the contraction of the Laplacian tensor 
A_: 

by a tensor <j, with the following definition: 

— (—A : a.a*p = F _1 [| (k,a.a*.k |$] = F _1 [| a* .k \ a ] (35) 



6 Levy case 



The second characteristic function of the increments SL of the (scalar) 
Levy forcing is the following: 

5K L (k, St) = StUki - D\k\ a (l - iP-^-u{k, a)] + o(5t) (36) 

|fc| 

where u(k,a) is defined by: 



71" Q/ 7f 

a ^ 1 : cu(k, a) = oj(a) = tan — ; a = 1 : u(k, a) = — log\k\ (37) 
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Considering an Ito-like forward integration of Eq.l, the increments 8L 
generates the following (first) characteristic function for the increments 
8X of the motion X(t): 

6Z x (k,6t\x-6x,t) = e tkm{x ' t] 5Z aL {k,5t\x,t) + +o(8t) (38) 
which yields the following elementary cumulant generating function SKx- 

SK x (k,8t\x,t) = 5t[ikm(x,t) +ik^a(x,t) (39) 

-D\k\ a (l - i^uik, a))a{x,t) a ] + o(5t) (40) 
|fc| 

and which is of the same type as Eq.17, with J = {1,q}. Therefore, as 
discussed in Sect. 5, we have fractional differentiations in the corresponding 
Eq.7, which will precisely correspond to Eq.4, and therefore establishes 
Prop. 1. 

Let us discuss briefly the regularity constraints that should be satisfied 
by the nonlinear function a((X(t),t) > and m((X(t), t). Obviously, they 
should be measurable. On the other hand, the uniqueness of the solution 
should require, as for the classical nonlinear Fokker-Planck equation (e;g. 
[18]), a Lipschitz condition for both a((X(t),t) and m((X(t),t). 



7 Extension to vector processes 

With but one important exception, the extension of the previous results 
to higher dimensions is rather straightforward. The starting point of this 
extension is the following nonlinear stochastic equation (2£(t) G R d ): 

dX{t) = m(X(t),t)dt + g{X(t),t).dL (41) 

where m and a are the natural vector, respectively tensor, extensions 
of the deterministic-like trend, respectively modulation of the random 
driving force. L is a d-dimensional Levy stable motion and, as discussed 
below, the expression of its characteristic function corresponds to the 
source of the difficulty in extending the scalar results to high dimensions. 
On the contrary, it is straightforward to check that Props. 2, 3 are valid 
in the d-dimensional case, with the following extensions (x £ R d ) for Eq. 
6: 

^(x,t\x ,t ) = J dy^(x-y\y,t)p(y,t\x Q ,t ) (42) 
and for Eq. 7 (n 6 J C N d , | n \= Y^ =1 rii): 

dp d 1 - 1 

-^{x,t\x ,t ) = 2^ q^I dx * d [ A n(2, 1 )pfe t \x , t )] (43) 

nej 1 2 " d 
the relation to the cumulants C n of the increments is now: 



9 



- (ni)!(n 2 )!..(n d )! - 

On the other hand, Eq. 41 yields the following extension to Eq.38: 

SZx(k, 6t\x, t) = e l ^^SZ^.iAk, 6t\x, t) (45) 
and therefore we have: 

5Kx{k, 6t\x, t) = i k.m(x, t) + 5Ki L (a*.k, 6t\x, t) + o(St) (46) 

Let us recall that a stable Levy vector in the classical sense [9, 15, 16] 
(see [17] for a discussion and a generalization) corresponds to the limit of 
a sum of jumps, with a power-law distribution, along random directions 
u 6 dB\, B\ being the unit ball, distributed according to a (positive) 
measure dS(u). The latter, which generalizes the scale parameter D of the 
scalar case, is the source of the difficulty since in general the probability 
distribution of a stable Levy vector depends on this measure, and therefore 
is a non parametric distribution. However, as discussed below, there is at 
least a trivial exception: the case of isotropic stable Levy vectors. 

Corresponding to our previous remarks, a (classical) stable Levy vector 
has the following (Fourier) cumulant generating function: 

Kl(Jc) = 5t[i{k,j) - / (ik,u) a dY;{u)] + o{St) (47) 

JuedBx 

which yields with the help of the Eq.46: 

^-Kx{k) = -div(m + g.j) - F^l / (ig(x,t).k,u) a dZ(u)] (48) 
The scalar case (Eq.36) corresponds to: 

0<p<l:/? = 2p-l, dT,(u) = Dcos( — )\p5 iu _ 1) + (l- P )5 (u+1) ] (49) 

For any dimension d, the second term on the right hand side of Eq.48 
corresponds to a fractional differentiation operator of order a. This op- 
erator can be slightly re-arranged. With the help of the odd dT,~ (u) and 
even dE + (u) parts of the measure dE(w), 

2 dE+(u) = dE(u) +dE(-u); 2 dE~ (u) = rfE(w) - dE(-u) (50) 
and the identity (0 being the Heaviside function): 

(ik) a = \k\ a [9(k)e 1 ^ + 6{-k)e- l£ ?] (51) 
one can write the extension of Eq.4 under the following form: 
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-g^p(x,t\x ,t ) = -div]m(x,t) + g(x,t).j)]p(x,t\x ,t )] 

-[< >e+ - < (S.-g*)-(-±: g-g*)** 1 >E-bfe%o>*o) (52) 

where the symmetric and antisymmetric operators are defined , similarly 
to Eq.35, in the following manner: 

- < i-A-g-g*)^ >E+ 

- < = g-g) 2 ^ >e- 

In general, each operator corresponds to a rather complex integra- 
tion (which is indicated by the symbol < . >e) of directional fractional 
Laplacians (Eq.33). However, the symmetric operator becomes simpler 
as soon as the even part d£ + of the measure d£ is isotropic. Indeed, the 
integration over directions yields only a prefactor D: 

< -(A '■ °L-£.*)^ >£+= D (-A : a-a*)% 

D= I d£+(u) | (u 1)M ) T (55) 

J uedBx 

and for a = 2 this corresponds to the classical term (A : <t.ct*) of the 
standard d-dimensional Fokker-Planck equation. If dS itself is rotation 
invariant, then the asymmetric operator vanishes, since <i£~ = 0. If 
furthermore, a is rotation invariant, i.e. a — al, then one obtains the 
following Fractional Fokker-Planck equation: 

d 

0^p(x,t\x ,t ) = -div\g.j_(x,t) + m(x,t)]p(x,t\x ,t ) (56) 
- D {( — A) a/2 ]a(x,t) a p(x,t\x Q ,to) (57) 

Therefore, as one might expect it, due to the rotation symmetries, this 
corresponds to a rather trivial extension of the standard gaussian case: a 
fractional power a of the d-dimensional Laplacian, as in the pure scalar 
case (Eq.4). Obviously, the integration performed in Eq.52 is also greatly 
simplified as soon as dE(w) is discrete, i.e. its support corresponds to a 
discrete set of directions u t . 

On the other hand, let us note that the framework of generalized stable 
Levy vectors [17], allows one to introduce a much stronger anisotropy than 
the the measure dS does if for classical stable Levy vectors. This there- 
fore diminishes the importance of the asymmetry of the latter. Indeed, 
the components of a generalized stable Levy vector do not have neces- 
sarily the same Levy stability index, the latter being generalized into a 



= / dS + (M)F _1 [| (g* {x,t).k,u) H (53) 

JuedB 1 

= I d^-(u)F- 1 [(-ig r (x,t).k,u) | (g*(x,t).k,u) \X5i] 
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second rank tensor. Similarly, the differential operators involved in the 
corresponding Fractional Fokker-Planck equation have no longer a unique 
order of differentiation. This is rather easy to check in case of a discrete 
measure dS(u) and we will explore elsewhere the general case. 

8 Conclusion 

We have derived a Fractional Fokker-Planck equation, i.e. a kinetic equa- 
tion which involves fractional derivatives, for the evolution of the proba- 
bility distribution of nonlinear stochastic differential equations driven by 
non-Gaussian Levy stable noises. We first established this equation in 
the scalar case, where it has a rather compact expression with the help of 
fractional powers of the Laplacian, and then discussed its extension to the 
vector case. This Fractional Fokker-Planck equation generalizes broadly 
previous results obtained for a linear Langevin-like equation with a Levy 
forcing, as well as the standard Fokker-Planck equation for a nonlinear 
Langevin equation with a Gaussian forcing. 
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